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ABSTRACT 

Standard models of radiation supported accretion disks generally assume that diffusive radiation 
flux is solely responsible for vertical heat transport. This requires that heat must be generated at 
a critical rate per unit volume if the disk is to be in hydrostatic and thermal equilibrium. This 
raises the question of how heat is generated and how energy is transported in MHD turbulence. 
By analysis of a number of radiation/MHD stratified shearing-box simulations, we show that the 
divergence of the diffusive radiation flux is indeed capped at the critical rate, but deep inside the 
disk, substantial vertical energy fiux is also carried by advection of radiation. Work done by radiation 
pressure is a significant part of the energy budget, and much of this work is dissipated later through 
damping by radiative diffusion. We show how this damping can be measured in the simulations, 
and identify its physical origins. Radiative damping accounts for as much as tens of percent of the 
total dissipation, and is the only realistic physical mechanism for dissipation of turbulence that can 
actually be resolved in numerical simulations of accretion disks. Buoyancy associated with dynamo- 
driven, highly magnetized, nearly-isobaric nonlinear slow magnetosonic fiuctuations is responsible for 
the radiation advection fiux, and also explains the persistent periodic magnetic upwelling seen at all 
values of the radiation to gas pressure ratio. The intimate connection between radiation advection 
and magnetic buoyancy is the first example we know of in astrophysics in which a dynamo has direct 
impact on the global energetics of a system. 

Subject headings: accretion, accretion disks — MHD — radiative transfer — turbulence 
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1. INTRODUCTION 

Much of the power radiated by the accretion fiow in 
luminous states of X-ray binaries and active galactic nu- 
clei necessarily originates from the vicinity of the central 
compact object where the gravity well is deepest. If the 
bolometric luminosity is anywhere close to Eddington, 
radiation pressure must dominate the thermal pressure 
of the accreting plasma in these regions, and understand- 
ing the physics of radiation dominated accretion is there- 
fore central to any explanation of how these sources work 
in t heir brightest states. Classical accretion disk the- 
ory (IShakura fc SunvaevlfTOTl INovikov fc ThorndflflTl 
models the fiow as geometrically thin and optically thick, 
and assumes that angular momentum is transported 
within the fiow by internal stresses due to turbulence. 
The average internal stress is postulated to scale with 
thermal, and therefore mostly radiation, pressure. Ac- 
cretion power is assumed to be locally dissipated as heat, 
which is then transported vertically outward by photon 
diffusion. 

This model has been questioned over the years on 
a number of grounds. First and foremost, if one 
takes its assumptions literally, then the resulting equi- 
librium structure is unstable to both thermal and in- 



fiow ( "viscous" ) ins t abilit ie s (iLightman fc Eardlevl 
oshil 



Shibazaki 
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19751: ISvunvaev fc Shakural 



1974; 



1975 



1976). However, these assumptions 



mav not be valid. IBisnovatvi-Kogan fc Blinnikovl 1 1977D 
point out that if the turbulent dissipation rate per unit 
mass is constant with height in the radiation pressure 
dominated regime, then the mass density would be in- 
dependent of height inside the photospheres, a situa- 
tion that would clearly be convectively unstableQ The 
argument is simple and worth repeating. Hydrostatic 
equilibrium in a radiation dominated disk implies that 
the vertial profile of radiation fiux F(z) with height z 
obeys kF/c = fl^z, where k is the opacity (which is 
very nearly constant as Thomson scattering dominates 
the fiux mean), c is the speed of light, il is the lo- 
cal angular velocity, and we have assumed Newtonian 

^ This assumption of constant dissipation per unit mass is im- 
pUcit in equation (2.22) of Sha|cura &_Sun tj.973), where they 
go beyond one-zone modeling and attempt a detailed treatment 
of the disk vertical structure. As a result, they conclude after 
their equation (2.26) that the density is independent of height. It 
has also been adopted by other authors even in recent years; e.g. 
appendix B of Begclman ( 2006) makes the same assumption in de- 
riving a simple model of convective transport in the presence of 
photon bubbles. 
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gravity for simplicity. If radiative diffusion dominates 
the heat transport, then equihbrium requires that the 
flux divergence dF/dz be equal to the dissipation rate 
per unit volume Q, giving Q = cVt} j n = Q* , which 
is constantly If the dissipation rate per unit mass is 
also assumed to be constant, then the density must be 
constant, implying convective instability. It should be 
noted, however, that vertical energy transport by con- 
vection under the conditions of constant dissipation per 
unit mass need not alter the thermal instability of the 
disk, as diffusive transport of energy could sti ll be dom- 
inant (jShakura. Sunvaev fc Zilitinkevichlll978| ). 

Energy may also be transported vertically in other 
ways. Rather than completely dissipating locally, mag- 
netic energy in the turbulence could be transported 
outward by buoyancy in the form of Poynting flux, 
perhaps to b e dissipated outside the disk ph otosphere 
in a corona (jGaleev. Rosner fc Vaianal I1979H . It has 
even been argued that magnetic buoyancy would limit 
the magnetic energy density in a radiation dominated 
plasma to be at most the gas pressure, thereby produc- 
ing an accretion stress that scales with the ga s pressure 
alone (iSakimoto fc CoronitilllQSlI : iStella fc Rosner, 1984 
iSakimoto &: Coronitilll989i r Such a stress would elimi- 
nate the thermal and inflow instabilities that plague the 
standard model. 

Unfortunately, this work was built on an incomplete 
foundation because it lacked an understanding of the 
physical nature of angular momentum transport in ac- 
cretion disks. One could therefore only guess the ver- 
tical profile of dissipation per unit mass, and pretend 
that the disk can be adequately modelled by a time- 
averaged, steady-state vertical structure. Moreover, ar- 
guments concerning the efficacy of magnetic buoyancy 
had to make assumptions about the rate of magnetic 
field generation, as well as the geometry of magnetic field 
lines. 

Since the discovery of the relevance of magne- 
torotation al (MR I) turb ul ence to accretion disks 
dlalbus fc Hawleyl [199TI : IHawlev fc Balb^ IT99TI : 
iBalbus fc Hawleyl 119981 ) . it has become possible to 
explore these ideas in detail with numerical simulations. 
The most relevant simulations thus far have used a 
shearing box that incorporates the vertical tidal gravi- 
tatio n al field of the cen tral object (Brandenburg et al. 
Il995t iStone et"ari 119961 ). The accretion stress in such 
simulations arises self-consistently from correlated 
magnetic and velocity fluctuations within the turbulence 
itself. Using such stratified shea ring boxes w i th an 
isothermal equation of state. Mi ller fc Stond ()200Clf ) 
found that the majority of the magnetic energy gener- 
ated in MRI turbulence was (numerically) dissipated 
locally. Nevertheless, a quarter of the magnetic energy 
generated was vertically transported outward by buoy- 
ancy, forming a strongly magnetized corona outside a 
weakly ma gnetized structure near the disk midplane. 
iMiller &rSt onc (2000) used an isothermal equation of 
state, however, and did not include the possibility of 
diffusive radiation transport outward along the vertical 
temperature gradients. 

This argument also led lShakura He Sunvae 3 {1971) to observe 
that 2cQ/{3k) is a characteristic value for the stress in radiation 
dominated disks. 



Inclusion of such transport would also introduce new 
dissipation physics, as compress ible waves should b e 
damped by radiative diffusion (|Agol fc KrolikI I1998D . 
This process is entirely analogous to the Silk damping of 
acoustic perturbations in the early universe (ISilH [19671 
11968.) . and to the radiative damping of stellar pulsation 
modes ()Coxlll980l ). It is particularly interesting because 
such dissipation can easily be resolved numerically in ra- 
diation MHD simulations. This contrasts sharply with 
our complete inability to resolve the microscopic scales 
on which viscous and resistive dissipation damps fluid 
and magnetic fluctuations. 

Shearing box simulations of MRI turbulence incor- 
porating thermodynamics and radiation transport have 
been possible for some time now. The first such simula- 
tions neglected vertical gravity and explored the prop- 
erties of the turbulence in the presence of radiation 
trans port, treated numerically using flux-limited di ffu- 
sion ((Turner. Stone fc Sand 120021: [Turner et al.ll2003[ ). If 
photon diffusion is rapid enough, and the magnetic pres- 
sure exceeds the pressure in the gas alone, then MRI tur- 
bulence can become extremely compressible with strong 
density fluctuations. These fluctuations are highly dissi- 
pative: net PdV work is done on the plasma over time, 
indicating an irreversible c onversion of mechanical en- 
ergy into internal energy ([Turner. Stone fc Sand [20021: 
[Turner et a l."2003') . This result conflrmed the suggestion 
of [Agol ~ Krolik (1998), although the speciflc character 
of the fluctuations being damped was not entirely clear. 

The flrst radiation MHD shearing box simulation of 
MRI turbulen ce with vertical gravity was published by 
[Turned ([2004[ ). The simulation was radiation pressure 
dominated, and exhibited no evidence of thermal insta- 
bility. Moreover, the time-averaged vertical entropy pro- 
flle was stable to hydrodynamic convection. On the other 
hand, the simulation was not fully energy-conserving, 
nor was it able to incorporate the photospheres within 
the simulation domain. Substantial mass loss occurred 
during the simulation and this was suggested as a pos- 
sible cause of the absence of any exponential thermal 
runaway. An attempt was also made to measure the ra- 
diative damping using the time-averaged PdV work, but 
the result was necessarily uncertain because this work 
can also be used to excite vertical mechanical motions 
when vertical gravity is present^ 

The technical issues with the 'Turner (2004") simulation 
were solved by Hirose, Krolik. fc Stone (2006). Radia- 
tion transfer was still treated by flux-limited diffusion, 
but a new diffusion solver permitted several improve- 
ments: The quasi-periodic radial boundary conditions 
appropriate to shearing boxes could be imposed prop- 
erly. Low densities could now be handled, allowing the 
photospheres to be incorporated within the simulation 
domain. Most importantly, a total energy scheme was 
implemented so that grid scale losses of magnetic and 
kinetic energy were fully captured as heat in the gas, 
and the code accurately conserves energy. While dis- 
sipation is therefore still numerical, it may nonetheless 
mimic local energy flow to microscopic dissipation scales 
at high wavenumbers in the turbulent cascade. Radiation 
and gas exchange momentum through Thomson scatter- 
ing and free-free opacity, and exchange energy through 
free-free absorption and emission. Energy exchange by 
Compton scattering was later incorporated into the code 
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in lHirose. Krolik. fc Blae"i (|2009l ). 

Using this improved code, we have succeeded in 
establishing long-lived therm al equilibria in which 
gas pressure dominates (jHirose. Krolik. fc St^cm^ 

f)6), ga s and r adiation pressure are comparable 
rolik. Hirose. fc Blaes ,20071). and radia tion pressure 
dominates ( Hirose. Krolik. fc BlaesI l2009f) . The radi- 
ation dominated state is thermally stable, even with 
no mass loss, and this is due to the fact that thermal 
pressure lags turbulent stress on time scales of order 
the thermal time (Hir ose. Krolik, & Blaes 2009). We 
have also established th at the stress and tot al pressure 
are correlated (see also lOhsuga et all I2009D . but only 
on time scales longer than the th ermal time, reviving 
the question of inflow i nstability ('L ightman fc EardlevI 
119741 iHirose. Blaes. fc Krolik 2009). Neither the gas 
density nor the dissipation rate per unit mass is constant 
with hei ght in the radiation dorninate d (or any other) 
regime (jHirose, Krolik. fc Blaei I20090 . and the time 
and horizontally-average d structures are stable to hy - 
drodynamic convection ([Krolik, Hirose, fc BlaesI [20071 ). 
However, t he outer la yers are generally magnetically 
supported fHirose. Krolik, & St one 200j), and exhibit 
Park er instability dynamics (jBlaes, Hirose, fc KrolikI 
[2001 . 

There remain a number of important questions about 
the thermodynamics of radiation dominated accretion 
disks. How does one properly calculate the contribu- 
tion of radiative damping to the overall dissipation when 
PdV work can also be used to excite vertical mechani- 
cal motions? What exactly are the compressive motions 
that are being dissipated by radiative diffusion? What 
happens if the local dissipation rate per unit volume ex- 
ceeds the radiation pressure dominated hydrostatic value 
of Q* — I n1 What controls the relative shares of ra- 
diative diffusion, advection, and Poynting flux in vertical 
energy transport, and how does this depend on radiation 
to gas pressure ratio and height? Can there be significant 
coherent energy flux in the form of, e.g., vertical acous- 
tic waves, and if so, how are these waves excited and 
how much energy do they transport and dissipate? Even 
if the time- and horizontally-averaged vertical structure 
is convectively stable, can local and transient buoyancy 
lead to significant energy transport? 

The goal of this paper is to answer these questions 
on the basis of detailed analysis of simulation data. In 
section [21 we provide a brief overview of the radiation 
dominated simulations we have analyzed. In section [3l 
we analyze global energetics, first deriving and discussing 
the total energy conservation equation in section lXTl and 
then turning to the first law of thermodynamics in sec- 
tion 13.21 There we show that work done by pressure is 
increasingly important at high radiation to gas pressure 
ratios, and this work is associated both with radiative 
damping and excitation of vertical mechanical motions. 
We identify two important classes of radiation pressure 
fluctuation in section [4l and show how radiative damping 
acts upon them in section[5l We then present detailed re- 
sults on the nature of vertical advective energy transport 
in section [6} In section [Tj we discuss how our findings 
give rise to a more dynamic view of the thermal physics 
of a radiation dominated disk. Finally, we summarize our 
conclusions in section [8] We provide some mathemati- 
cal background on trapped vertical modes that modulate 



the mechanical work and advective energy transport in 
an appendix. 

2. SIMULATIONS 

The radiation MHD equations solved in our simu- 
lations are discussed by Hirose, Krolik, fc Stond ([20061 ) 
and fHirose, Krolik, fc Blaesi ([2009 ), but we list them here 
again as they are the basis for most of the equations we 
derive elsewhere in the paper. 



dt 



V ■ (pv) = 



(1) 



d 1 
— {pv) + V- (/9vv) = -Vp - V - q + —(V X B) X B 
at 47r 

+ + fsB (2) 

c 



de 
dt 



V • (ev) : 



-pV ■ V - q : Vv - (aT^ - E)cK^p 



-cEkcsP- 



^ + V • (£;v) = -P : Vv + {aT^ - E) 



Q 



CKsP 



(3) 



4fcB(T - Trad) „ „ / 

-cEk^sP - V • F (4) 



'dt 



= V X (v X B) 



cA 
R^p 



VE 



(5) 
(6) 



Here p is the density, v is the fluid velocity, B is the 
magnetic field, p is the pressure in the gas, q is a diag- 
onal tensor associate d with the artif icial bulk viscosity 
adopted by the code ([Stone fc Norma n 1992), e = 3p/2 
is the internal energy density in the gas, T is the gas 
temperature, E is the radiation energy density, P is the 
radiation pressure tensor. Trad = {E/a)^/^ is the effective 
temperature of the radiation, a is the radiation density 
constant, F is the radiation flux, is Boltzmann's con- 
stant, TOg is the electron mass, Kqs is the electron scatter- 
ing opacity, = (p, e) is the Planck mean free- free 
opacity, kF^ = K^(p, e) is the Rosseland mean opacity in- 
cluding electron scattering and free- free contribution^, A 
is a flux limiter equal to 1/3 in the optically thick limit, 

and Q is the dissipation rate per unit volume required 
to maintain total energy conservation due to grid-scale 
losses of magnetic and kinetic energy. More details on 
how many of these quantities are actually computed in 
the code can be found in|Hirpse. Kr olik. & Stone (200(f) 
and IHirose. Krolik. fc BlaesI ([20091 ). 

The vector fss represents the gravitational and inertial 
forces in the local shearing box frame, which is rotating 
at fixed angular velocity 51 with respect to an inertial 
reference frame: 



fsB = —ipflz XV- pV0, 



(7) 



^ Electron scattering dominates the Rosseland mean opacity in 
all the simulations considered in this paper, so that ft'' = Kes quite 
accurately. 
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where 4> = —31)^x^/2 + f2^z^/2 is the effective gravita- 
tional potential, and x, y, z are Cartesian spatial coordi- 
nates along the radial, azimuthal, and vertical directions, 
respectively. 

Table [1] summarizes most of the numerical parameters 
of the six simulations we analyze in this paper, includ- 
ing the surface mass density S and scale height H of the 
initial condition. Each simulation had the same angular 
velocity Vl — 190 rad s~^, corresponding to a radius of 
iQGM/c^ around a 6.62 Mq Schwarzschild black hole. 
The grids all had 48 zones in the x (radial) direction, 96 
zones in the y (azimuthal) direction, and 896 zones in the 
z (vertical) direction. The simulations were all initialized 
with a weak magnetic field in a twisted azimuthal flux 
tube geometry in the midplane regions of the domain. 
Magnetorotational turbulence is well-established by ten 
orbits into the simulation, and we ran each simulation a 
further 250 to 600 orbits. This is much longer than the 
thermal time of the simulations, the average of which 
ranges from 13 orbits for 1112a to 24 orbits for 0519b, 
as shown in Table [2] All the simulations reach approx- 
imate thermal equilibria with continued long time scale 
fluctuations. Radiation pressure dominates the ther- 
mal pressure in all these simulations, and the time and 
box-averaged ratios of radiation to gas pressure are also 
listed in Tabled] M ore information about these sim ula- 
tions can be found in 'Hirose . Krolik. fc BlaesI (|2009D and 
iHirose. Blaes. fc Kr olik (200^^ 

3. ENERGETICS AND THERMODYNAMICS 

The shearing box equations of motion as implemented 
in our code conserve total energy to high accuracy 
(jHirose, Krolik. fc Stonell2"006D . Energy originates in the 
work done by turbulent stresses on the shearing radial 
walls of the box and ultimately escapes from the top and 
bottom of the box, largely in the form of photons. Long- 
term equilibrium requires that nearly all the work done 
by the walls be dissipated and that all the heat gener- 
ated be vented. In this section we analyze this energy 
flow and dissipation in detail as a function of height. 

3.1. Total Energy Conservation 

Equations (P)-® and ([7]) together imply a total energy 
conservation equation 

d 

"q^ (*^incch ~(" "^thcrm) ^ ' (Fmcch ~l~ Ftiicrm) 

= v (^V-P + ^F^ . (8) 

The grid scale numerical heating Q has been lost, as our 
total energy numerical scheme is designed to capture grid 
scale losses of magnetic and kinetic energy and convert 
them into heat, which is promptly used to create photons. 
There are also numerical energy sources and sinks that 
should be present on the right hand side of this equation 
due to density and i nternal energy floors and a ve locity 
cap (see appendix of lHirose. Krolik. fc Stonell2006[ ). but 
these are negligible and we ignore them here. 

The energy density in equation ([8]) has two parts. The 
first is mechanical, being the sum of kinetic, effective 
gravitational potential, and magnetic energy densities: 



.'mech 



Stt' 



(9) 



The second is thermal, being the sum of the gas and 
radiation internal energy densities: 

£thcrm = e + E. (10) 

The energy flux vector also has two similar pieces. The 
mechanical energy flux is the sum of kinetic energy fiux, 
effective gravitational potential energy flux, flux of work 
done by artificial viscosity and by gas and radiation pres- 
sures, and Poynting flux: 

1 c 
Fmoch = -pv^v + p0v + q-v+pv+P-v+— ExB, (11) 
Z 47r 

where E = — v x B/c is the electric field in ideal MHD. 
The thermal energy flux is the sum of gas and radiation 
internal energy advection and heat transport by radiative 
diffusion. 



'the 



E)w- 



(12) 



The right hand side of equation ^ is an artificial set 
of energy source and sink terms that result from our use 
of flux- limited diffusion to handle radiation transport. 
These terms exactly cancel in the optically thick limit. 
They would also cancel to lowest order in v/c if we were 
using the full radiation moment um equation (eq. [9] of 
IStone. Mihalas fc Normanlll992D instead of the diffusion 
equation ([6]). These terms contribute negligibly to the 
overall energy balance in our simulations, and we ignore 
them from now on. 

Assuming an equilibrium has been established over 
long time scales, the energy conservation equation ([8]) 
can be rewritten to show that the local divergence of the 
time and horizontal average of the vertical energy flux is 
given by the local rate at which work is being done on 
the fluid by the shearing walls at that height z, i.e. 

d 3 

^ ((i^mcch,.) + (Fthorm,.)) = 2^ i^^vi^)) ' (1^) 

where angle brackets denote horizontal and time- 
averages at a particular height z. For example. 



1 



At rLy/2 

dt I dx d?/i^mcch,^(14) 

L^/2 



where At is the simulation duration minus the first 10 
orbits when the MRI was still in its growth phase. All 
time-averages presented in this paper use this time inter- 
val. 

The quantity {Txy{z)) is the time-averaged xy compo- 
nent of the Reynolds and Maxwell stresses at height z, 
plus a negligibly small contribution from radiation vis- 
cosity. 



{Txy{z))--- 



1 



AtL 



At 



rLy/2 



V Jo 

BxBy 

47r 



dt 



P. 



-Ly/2 



dy\ pVxSvy 



(15) 



Here Sv,, 



3Qx/2 is the perturbation of azimuthal 



velocity from the average shear flow, and the y-integral 
is d one either on the inner or outer radial wall of the 
box (|Hawlev. Gammie. fc Balbu£lll995l) . When vertically 
integrated, the relative contributions of Maxwell and 
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TABLE 1 

Simulation Parameters 



Simulation S 


H 


Duration 


Box Dimensions 


(g cm-"^) 


(cm) 


(orbits) 


(L^/H X Ly/H X L,/H) 



5.43 X 10'' 5.83 x 10^ 264 0.3375 x 1.35 x 6.3 

7.48 X lO* 4.37 x 10'^ 403 0.3375 x 1.35 x 6.3 

1.06 x 10^ 1.46 X 10"^ 610 0.45 x 1.8 x 8.4 

1.06 X 10^ 1.46 X 10"^ 611 0.45 x 1.8 x 8.4 

1.24 X 10^ 1.17x10'^ 603 0.54 x 2.16 x 10.08 

1.52 X 10^ 7.28 X 10^ 426 0.6 x 2.4 x 11.2 



TABLE 2 

Summary of Time-Averaged Physical Properties of 
Simulations 



Simulation Hp/H^ f thermal /orbits^ {^^^Y (j^ 



0211b 


0.55 


16 ±3.2 


62± 16 


2.3 ±0.56 


0519b 


0.81 


24 ±5.2 


70 ±21 


2.6 ±1.3 


1112a 


0.71 


13 ± 2.6 


6.6 ±2.4 


1.9 ±0.54 


1126b 


0.88 


16 ±2.9 


10 ± 2.7 


2.1 ±0.81 


0520a 


1.3 


22 ± 5.3 


14 ±4.7 


2.2 ±0.79 


0320a 


1.5 


21 ± 3.8 


6.6± 1.6 


1.7 ±0.33 



^ The thermal pressure scale height, defined as Hp = 
f^oo ^therm(^)'^2/[2Ptherm(0)], where PthermCz) is the time-avcraged 
vertical profile of gas plus radiation pressure, in units of the fiducial 
scale height H used in the simulation grid. 

^ The time-average of the thermal time, defined as the ratio of in- 
stantaneous thermal energy content in the simulation domain to the 
instantaneous horizontally averaged radiative flux emerging from the 
top and bottom faces. Errors indicate one standard deviation in the 
time average. 

The time average of the ratio of box-averaged radiation pressure 
E/3 to box-averaged gas pressure p. Errors indicate one standard 
deviation in the fluctuations of the ratio about the time average. 

The time average of the ratio of maximum to minimum density at 
the midplane z = 0. Errors indicate one standard deviation in the 
fluctuations of the ratio about the time average. 

Reynolds stresses are remarkably constant from simu- 
lation to simulation: they are 85% and 15% respectively, 
to within 1% for all six simulations. The radiation vis- 
cosity contributes at most 10""' of the total vertically 
integrated stress. 

Figure [1] shows the time-averaged vertical profiles of 
the dominant contributions to the vertical thermal and 
mechanical energy fluxes in simulations 1112a (one of 
the two lowest radiation to gas pressure ratio simula- 
tions we consider in this paper) and 0519b (the highest 
radiation to gas pressure ratio) . As in the case of gas 
press ure dominated simulations (jHirose. Krolik. fc Stoi^ 
|2006( ) and simulations with comparable gas and radiation 
pressure (Krolik, Hirose. fc Blaes 2007.), radiative diffu- 
sion is the dominant process of vertical energy trans- 
port in simulation 1112a. However, we now see that at 
the highest levels of radiation pressure support simulated 
thus far (0519b), radiation advection is just as important 
in the midplane regions. 

Figure [2] depicts the time-averaged vertical profiles of 
the various contributions to the divergences of the ver- 
tical thermal and mechanical energy fluxes, and com- 
pares their sum to the vertical profile of stress times 
rate of strain. Equation (jl3|) is accurately satisfied by 
all the simulations. The radiation diffusion flux diver- 



0211b 
0519b 
1112a 
1126b 
0520a 
0320a 



1.0 




-4 -2 2 4 




-3-2-10 1 2 3 

z/H 



Fig. 1. — Vertical profiles of the most important contributors to 
the horizontally and time-averaged vertical energy fiux in simu- 
lations 1112a (top) and 0519b (bottom). The different curves are 
diffusive radiation energy flux (red) , advected radiation energy flux 
(dashed red), flux of radiation pressure work (gray), Poynting flux 
(blue), and advected gas internal energy flux (green). The radi- 
ation pressure work and Poynting flux curves nearly coincide in 
simulation 1112a. 

gence approximately matches the cVi^ / Kcs value required 
by hydrostatic equilibrium, presumably because depar- 
tures from this equilibrium would result in very fast read- 
justments on the dynamical time scale. All the remain- 
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ing flux divergence components are all positive in the 
midplane regions and negative further out. The energy 
injected by the stresses on the walls in the midplane re- 
gions is significantly larger than can be carried away by 
radiative diffusion and still maintain vertical hydrostatic 
equilibrium. This excess energy is therefore transported 
outward by the other forms of energy flux, among which 
radiation advection is the most important. This process 
is generally completely ignored in standard accretion disk 
models. 

3.2. The First Law of Thermodynamics: PdV Work 
and Radiative Damping 

By combining equations ([3]) and Q, energy conserva- 
tion may also be described in terms of the first law of 
thermodynamics, 

+ V . Fthcon = g - pV • V - P : Vv. (16) 

ot 

Here Q = Q — q : Vv is the local dissipation rate per unit 
volume, being a sum of grid-scale conversions of magnetic 
and kinetic energy into heat as well as energy dissipated 
by artificial viscosity. (The artificial viscosity dissipation 
is approximately 3/4 of the grid-scale loss rate of kinetic 
energy in all of our simulations, and has a similar time- 
and horizontally-averaged vertical profile.) Again assum- 
ing an equilibrium has been established over long time 
scales, equation ([TE)) implies that 

^ (Fthcrm,.) = (Q - pV • V - P : Vv) . (17) 
dz 

In other words, the local divergence of the average ver- 
tical thermal energy fiux is equal to the local dissipation 
per unit volume plus the rate at which gas and radia- 
tion pressure work is being done on the plasma per unit 
volume. 

Equation ([T7]) appears to be familiar, but its form hides 
some subtleties peculiar to radiation dominated disks. 
These new effects emerge when comparing the two sides 
of this equation. Figure [3] depicts the vertical profiles of 
the advection and radiative diffusion contributions to the 
divergence of the vertical thermal flux, as well as the nu- 
merical dissipation profile, for all six of our simulations. 
Because our numerical dissipation Q does not include 
damping by radiative diffusion, one might have expected 
it to be smaller than the divergence of the thermal fiux. 
In fact the opposite is true, with the largest discrepancy 
occurring in the most radiation pressure dominated sim- 
ulation 0519b. Moreover, the absence of an explicit place 
for radiative damping suggests that something is missing 
because that process should contribute to the total dis- 
sipation rate. 

Figure|3]shows the same vertical profiles of thermal flux 
divergence, but a dissipation rate corrected by the PdV 
work. With this adjustment, the time-averaged first law 
of thermodynamics ([TT]) is accurately satisfied in each of 
the simulations. From this fact, we reach two impor- 
tant conclusions. First, there is a significant conceptual 
flaw in classical time-steady accretion disk models. These 
models assume that the divergence of the diffusive radi- 
ation flux completely defines the left hand side of the 
first law of thermodynamics, while the local dissipation 
rate is the only contribution to its right hand side. But 



we have just seen that in radiation dominated disks both 
of these simplifications are wrong: radiation advection 
must be included with radiative diffusion, and the PdV 
work terms are important. Second, radiative damping is 
actually included in the PdV work terms. In our scheme, 
the dissipation associated with it is conveyed by the dif- 
fusion equation, but the energy it dissipates fiows into 
the gas via the pressure work terms in the first law of 
thermodynamics. 

It is perhaps helpful to make more explicit the actual 
dissipation associated with radiative diffusion. The first 
law of thermodynamics (|16p can be combined with the 
radiative diffusion equation (jS]) to derive an equation for 
the evolution of the entropy per unit mass s of the gas 
plus radiation mixture. Restricting consideration to op- 
tically thick regions for simplicity, this equation is 

f ds \ /F\ AacT , ,9 Q , , 

The last term on the left hand side is the divergence 
of the entropy flux due to radiative diffusion. The flrst 
and second terms on the right hand side are sources of 
entropy (dissipation) due to radiative diffusion and grid 
scale losses plus artiflcial viscosity, respectively. Assum- 
ing an equilibrium has been established on long time 
scales, this equation becomes 

(19) 

Advective and diffusive transport of entropy is therefore 
balanced by dissipation. The dissipation associated with 
radiative damping of fluctuations in the turbulence is in 
the first term on the right hand side, but this term also 
includes entropy generation due merely to the photon 
diffusion down the background average vertical temper- 
ature gradient. An analogous situation holds in stars 
with radiative envelopes. Under static conditions, the 
outward luminosity at every radius in the envelope is 
constant. Because the temperature at the base of the 
envelope is higher than the temperature of the photo- 
sphere, the entropy leaving the photosphere exceeds the 
entropy entering the base of the envelope. The source of 
this entropy increase is the dissipative nature of photon 
diffusion itself, and is not associated with dissipative re- 
lease of energy that must then be transported away to 
establish thermal equilibrium. 

Because of this non-energy releasing dissipation, we 
cannot use equation (fT9|) to calculate the true dissipa- 
tive heating due to radiative damping of fluctuations. 
Instead, we must try and use the pressure work terms in 
the flrst law of thermodynamics P7|) . As wc noted above, 
it is through these terms that energy flows into the gas 
from turbulent fluctuations, and ultimately dissipates by 
radiative diffusion. 

However, somewhat surprisingly, the time-averaged 
pressure work corrections — < pV ■ v > — < P : 
Vv > are typically negative., in contrast to the case 
of shearing boxes without v ertical gravity stud ied by 
[Turner. Stone fc Sand (|200l : [Turner et all (|200l . This 
sign change indicates that the plasma does net work 
through expansion even though these corrections must 
include radiative damping, which would be positive (it 
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Fig. 2. — Divergences of the horizontally-averaged and time-averaged vertical thermal and mechanical energy fluxes. The different colored 
curves show the divergences of various vertical fluxes: advected gas internal energy flux (green), advected radiation energy flux (dashed 
red), diffusive radiation energy flux (solid red), flux of radiation pressure work (gray), and Poynting flux (blue). The total of these flux 
divergences is shown as the solid black curve, which matches very well the proflle of stress times rate of strain (dashed black curve), in 
agreement with equation JTHJ . (We have neglected the flux of gas pressure work, which being 2/3 of the advected gas internal energy flux 
is negligible.) The horizontal dotted line indicates the fiducial dissipation cC^/kos expected from hydrostatic equilibrium in the radiation 
pressure dominated limit. 



is dissipative) . The fact that the overall pressure work 
corrections are negative means that the gas loses more 
thermal energy by driving expansion away from the mid- 
plane than it gains by radiative damping. 

Fortunately, these two contributions to the pres- 
sure work corrections (pumping expansion and radia- 
tive damping) can be separated analytically. Consid- 
ering only the optically thick regions, write the ther- 
mal pressure Pthcrm = p + E/3 and fluid velocity as 
-Ptherm,av + <5Ptherm and Vav + Sv , respectively. Here the 
"av" subscript denotes horizontal average and ^Ptficrm 
and (5v therefore have zero horizontal average by defi- 
nition. The divergence of also has zero horizontal 
average, as one can show by integrating by parts and 
using the shearing box boundary conditions. Hence the 
horizontal average of the thermal pressure work done on 
the plasma under optically thick conditions is 

I /•Li,/2 rLy/2 

/ dx dy {-pV • V - P : Vv) = 

LxLy J-Ly/2 

d 

-Pthcrin.av ^av,z(^5 ^) 

oz 

+ —— dx dy(-<5PthermV-(5v).(20) 

^x^^y J-L^/2 J-Ly/2 

After time-averaging, the first term on the right hand 
side represents pumping of vertical mechanical motions. 



This term will be negative if the plasma undergoes ver- 
tical expansion in a horizontally-averaged sense, i.e. the 
plasma will do net work. We will henceforth call this 
term the mechanical pumping term. 

The second term on the right hand side of equation 
(pn|) is the radiative damping. For a sinusoidal adiabatic 
fluctuation, this term would vanish identically in the time 
average, because SPtham and V • 6v are 90 degrees out 
of phase. Radiative diffusion removes this cancellation 
by making iJPtherm lead V • (5v by a little more than 
90 degrees. For example, if the plasma is locally com- 
pressed nearly adiabatically, the temperature rises and 
radiative diffusion to the cooler surrounding regions in- 
creases. This diffusion is fastest as the point of maximum 
compression is approached, while at the same time the 
rate of compression is slowing down. The temperature 
therefore starts to drop just before the point of maxi- 
mum compression, i.e. the maximum in temperature is 
reached before the point of maximum compression. The 
same holds on the expansion cycle, where the minimum 
in temperature is reached before the point of maximum 
expansion. Hence there is a greater than 90 degree lead 
in ^Pthcrm relative to V • (JvO This phase offset results 

* Very short wavelength fluctuations are approximately isother- 
mal, not adiabatic, because radiative diffusion is so rapid on short 
length scales. Finite (as opposed to infinitely rapid) radiative dif- 
fusion in this limit also produces the same greater than 90 degree 
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Fig. 3. — Divergences of the horizontally-averaged and time-averaged vertical thermal energy fluxes. The green curve shows the divergence 
of the advected gas internal energy flux, the dashed red curve shows the divergence of the advected radiation energy flux, the solid red 
curve shows the divergence of the diffusive radiation energy flux, and the solid black curve shows the total of these three curves, i.e. 
(di^thcrm,z/'^^)- The dashed black curve shows the horizontal and time-averaged dissipation rate per unit volume < Q >. This exceeds 
the divergence of the thermal flux, the most notable discrepancy being in simulation 0519b, which has the highest average radiation to gas 
pressure ratio. The horizontal dotted line indicates the fiducial dissipation cH^/kos expected from hydrostatic equilibrium in the radiation 
pressure dominated limit. 



in a net positive time average of — (5-PthormV ■ ^v. The 
plasma has net positive work done on it, and that excess 
work is dissipated by radiative diffusion. In principle, 
the mechanical pumping term of equation ()20p may still 
include a smah positive contribution from the radiative 
damping of purely vertical acoustic waves, but we show 
in section [01 below that this contribution is neghgible. 

Figures [S] and [B] depict the time-averaged vertical pro- 
files of stress times rate of strain and the different contri- 
butions to dissipation and mechanical work for simula- 
tions 1112a and 0519b. We did not save high time reso- 
lution data over the entire simulation duration to enable 
us to directly compute the last term on the right hand 
side of equation ()20p . so the radiative damping profiles 
in these figures were computed from subtracting the ver- 
tical profile of the first term from the vertical profile of 
the left hand side. We also neglected gas pressure work 
here as it is very small in these simulations. 

Two features are worth noting about these profiles. 
First, even after time-averaging, there remain spatial 

lead between thermal pressure and V • (5v. In a compression phase 
in this case, the work being done on the plasma causes it to be 
a little hotter than it would be if it were isothermal, resulting in 
a slightly greater thermal pressure. But this excess pressure must 
then drop as the point of maximum compression is reached because 
radiative diffusion is most rapid there, returning the pressure to the 
isothermal value. Hence the maximum in thermal pressure leads 
the point of maximum compression. 



fluctuations in the total mechanical work proflle as well 
as the profile of mechanical pumping. Those fiuctuations 
are completely absent in the radiative damping profile, 
which is as smooth as the other (magnetic and kinetic) 
dissipation profiles. The radiative damping profile is also 
very similar in shape to these other dissipation profiles. 
Second, when compared to the numerical dissipation, the 
radiative damping and mechanical pumping, as well as 
the net total pressure work, are clearly relatively more 
important in 0519b (the simulation with the highest ra- 
diation to gas pressure ratio) than in 1112a. As much as 
22 percent of the work done by the shearing walls ends 
up being dissipated by radiative diffusion in simulation 
O519b0 

Figure [7] illustrates this trend of increasing relative im- 
portance of the pressure work corrections with growing 
radiation to gas pressure ratio. A larger rate of radiative 
damping presumably requires larger density fluctuations; 

^ The radiative damping percentages of the work done by the 
shearing walls are 12 and 15 percent for simulations 1112a and 
1126b, respectively. These are much higher than the 1.3 and 0.7 
percent values that we stated in 'Hirosc, Krolik, & Blacs (2003) (see 
end of section 3 in that paper). Those previous numbers came from 
integrating up the total pressure work near the midplane where it 
is (barely) positive (see the solid gray curve of FigurefS)!, indicating 
net damping. Our new analysis, which cleanly separates radiative 
damping from work associated with vertical expansion, shows that 
the true radiative damping is much larger in these simulations and 
peaks off the midplane (upper gray dashed curve of Figure [Sj. 
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Fig. 4. — Same as Figure |3] except now the dashed black curve shows the horizontal and time-averaged dissipation plus the rate at which 
pressure work is being done on the plasma per unit volume, i.e. < Q — P : Vv >. (The gas pressure work — pV ■ v is negligible, and has 
been neglected.) This agrees very well with the divergence of the thermal flux, in agreement with equation l|17p . The horizontal dotted 
line indicates the fiducial dissipation cQ? / k,qs expected from hydrostatic equilibrium in the radiation pressure dominated limit. 
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Fig. 5. — Time and horizontally-averaged vertical profiles of 
stress times rate of strain 3f2T^^/2 (dashed), grid scale mag- 
netic energy dissipation (blue), grid scale kinetic energy dissipa- 
tion plus artificial viscosity dissipation (green), and — P : Vv work 
(grey), for simulation 1112a. The total of these last three, i.e. 
Q — P : Vv, is shown as the solid black curve and matches the 
time and horizontally-averaged profile of thermal energy fiux di- 
vergence. The lower gray dashed curve shows the time-averaged 
profile of —{E!iv/3)dvz,!i,v/dz, i.e. minus the horizontally averaged 
radiation pressure times the vertical derivative of the horizontally 
averaged vertical velocity. This represents the spatial profile of 
work done to pump vertical mechanical motions. The difference 
between this and the total pressure work profile is given by the 
upper dashed gray curve, which represents the radiative damping 
contribution to the dissipation. The horizontal dotted line again 
indicates the fiducial dissipation cH^/kos for a radiation pressure 
dominated hydrostatic equilibrium. Vertical dotted lines indicate 
one pressure scale height Hp away from the midplane. 
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Fig. 6. — Same as Figure [5l except for simulation 0519b. 

to test this supposition, we have measured the time aver- 
age of the ratio of maximum to minimum density at the 
midplane in each of the simulations, and the results are 
listed in Table [2] Density fluctuations in the midplane 
regions do indeed become slightly larger with increasing 
radiation to gas pressure. Figure [8] shows the fractional 
pressure work as a function of the time-averaged density 
contrast at the midplane. The large horizontal error bars 
reflect the large variations in the density contrast, but the 
radiative damping points of this figure are clearly consis- 
tent wit h the trend observe d in non-stratified shearing 
boxes by [Turner et al.l (|2003[ ) (see top panel of their Fig- 
ure 7). 
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Fig. 7. — Time-averaged ratios of various box-averaged contri- 
butions to the dissipation and pressure work to the vertically in- 
tegrated stress times rate of strain, as a function of the average 
radiation to gas pressure ratio of each simulation. The upper set 
of points (crosses) shows the fractional contribution of the numer- 
ical dissipation Q (grid scale losses of magnetic and kinetic energy 
as well as artificial viscosity). The middle set of points (diamonds) 
shows the radiative damping, and the bottom set of points (trian- 
gles) shows the radiation pressure work associated with pumping 
of vertical mechanical motions. Error bars indicate one standard 
deviation in the time-averages. 
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Fig. 8. — Same as Figure [T] only now as a function of the time- 
averaged density contrast at the midplane of each simulation. 

There is considerable structure in the temporal behav- 
ior of the pressure work terms. A ten orbit segment of 
the time-dependence of the vertical integral of the stress 
times rate of strain, the dissipation terms, and the pres- 
sure work terms for simulation 0519b is shown in Fig- 
ure [SI The mechanical pumping work (lower gray dashed 
curve) shows clear oscillatory behavior on time scales of 
order the orbital period. This is completely absent in the 
radiative damping (upper gray dashed curve), which in- 
stead clearly exhibits much higher frequency variability. 
Both also exhibit much longer time scale variation. 

Some aspects of this behavior can be immediately 
understood by Fourier transforming the time depen- 
dence and plotting the temporal power spectrum of the 
vertically integrated mechanical pumping and radiative 
damping terms. The result is shown in Figure 1101 The 
large oscillations seen in the mechanical pumping in Fig- 
ure [9] are reflected in a series of discrete sharp peaks. 
These peaks represent standing vertical acoustic waves 
that are trapped in the box. Despite their prominence, 
we show below in section 16.11 that they actually con- 
tribute negligibly to the energetics of the disk. The ra- 
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Fig. 9. — Time dependence over a ten orbit period of the verti- 
cally integrated stress times rate of strain and various contributions 
to the dissipation and work shown in Figure|5]for simulation 0519b. 
The colors and line styles correspond to the same quantities as in 
Figures \5\ and [S] 




0.01 0.10 1.00 10.00 

f /(orbital period) ' 

Fig. 10. — Vertically-integrated power spectra of the — P : Vv 
work (black), the mechanical pumping portion of this work (green), 
and the radiative damping portion of this work (red), for simulation 
0519b. 

diative damping power spectrum exhibits much higher 
frequency power, including a set of high frequency, ex- 
tremely sharp spikes. This is consistent with the high 
frequency variability seen in the time domain of Figure |9l 
To better understand this variability, we must first ex- 
amine the types of fluctuation that can produce radiative 
damping, to which we now turn. 

4. VARIETIES OF RADIATION PRESSURE 
FLUCTUATIONS 

It is radiation pressure (i.e. temperature) fluctuations 
that give rise to photon diffusion and therefore radiative 
damping. In this section we discuss two distinct types of 
such fluctuations that clearly play a role in our simula- 
tions. We will also see in section [HI that understanding 
both radiation advection and energy transport by Foynt- 
ing flux will likewise be aided by a prior understanding 
of radiation pressure fluctuations. 

In classical MHD theory (in which radiation pressure is 
negligible), linear compressible waves are classified into 
"fast" and "slow" modes. One way to understand quali- 
tatively this distinction is to note that the magnetic and 
gas pressure perturbations are exactly in phase in the 
former mode and exactly out of phase in the latter; it 
is the partial cancellation of the pressure perturbation 
in the slow mode that causes its low propagation speed. 
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Our simulations do not exactly correspond to this cate- 
gorization in two ways: the fluctuations are often nonlin- 
ear; and radiation pressure is both significant and only 
imperfectly coupled to the fluctuations (that is what ra- 
diative damping is all about, of course). Nonetheless, 
this conceptual division remains useful because large to- 
tal (gas plus radiation plus magnetic) pressure fluctua- 
tions at a given length scale tend to have much higher 
frequency than those with small total pressure fluctu- 
ations. The large total pressure fluctuations therefore 
propagate on an approximately stationary background 
set by the slowly evolving small total pressure fluctua- 
tions. We will call the former "acoustic waves" and the 
latter "isobaric waves" in order to stress this distinction. 

4.1. Acoustic waves 

lAgol &: KrolikI (|1998[ ) suggested that acoustic waves 
(i.e., fast magnetosonic waves) would be the dominant 
contributor to radiative damping in radiation dominated 
accretion disks. We have already seen that standing 
vertical acoustic waves are excited in the box. In ad- 
dition, strong nonaxisymmetric acoustic waves are per- 
vasive in all our simulations. These waves are almost cer- 
tainly stochastically excited by th e MRI turbulence itself 
(jHeinema nn fc Papaloizoull2009allH ) . In these waves, gas 
density, radiation pressure, and magnetic pressure fluctu- 
ations all oscillate in phase. Figure [TT] depicts a snapshot 
of various quantities at the z — midplane of simulation 
0519b at 200.3 orbits. A wave pattern is clearly evident 
in most of the fluid quantities shown, though it is clean- 
est in the total pressure (upper left). Figure [T^ depicts 
a vertical slice of the midplane regions at the same time. 
Particularly in the total pressure (upper left), the wave is 
once again evident in the vertical yellow and blue stripes 
in the midplane region \z\ <^ 0.2 — 0.3H and primarily 
propagates in the horizontal direction. 

To better understand these nonaxisymmetric waves, it 
is helpful to project the spatial variation of the z = 
midplane fluid quantities onto the natural set of basis 
vectors exp[i{kx{t)x + kyy)] of the shearing box, where 

27rn 3 

kccit) = -j-^ + -nkyt, (21) 



and rix a nd Uy are integer qua ntum num- 
bers (iHawlev. Gainmie. fc Balbui 119951 : 

iHeinemann fc Papaloizoul l2009b[ ). The largest az- 
imuthal wavelength {\ny\ = 1) waves always have the 
most power. Figure [T3] shows the real part of the total 
pressure Fourier amplitude as a function of time for 
n„ = 1 and various valu e s of ri x- In agreement with 
IHeinemann fc Papaloizoul (|2009bD . waves with different 
values of the quantum number reach high amplitude 
at distinct times (in declining order) when they swing 
from leading to trailing, resulting in a series of wave 
pulses. The separation in time between these pulses 
is determined entirely by the time interval between 
successive epochs at which the shearing radial bound- 
aries become exactly periodic in the radial direction. 
This shear time of the box is STs = 2Ly/{3flLx), i.e. 
37r/4 ~ 2.36 inverse orbital periods. We purposely 



chose the 200.3 orbits epoch for Figures [Tl] and [T^ as 
it corresponds to the time of peak amplitude for one of 
these wave pulses {rix — 9, Uy — 1). 

4.2. Isobaric waves 

The nonaxisymmetric acoustic wave pattern is evident 
in the spatial distribution of density, shown in the top 
right panel of Figure [Tl] but this pattern is markedly per- 
turbed by shorter length scale fluctuations. Among the 
most prominant are rarefied regions, e.g. near {x — 0.05, 
y — —0.5), that are correlated with regions of low ra- 
diation pressure, low gas pressure, and high magnetic 
pressure (bottom left to right panels of Figure [TTl re- 
spectively). 

These are a second kind of radiation pressure fluctua- 
tion, one in which the magnetic pressure oscillates with 
very nearly the same amplitude as the sum of gas and 
radiation pressure, but with opposite phase. As a result, 
the total pressure hardly changes. This near-cancellation 
of pressure fluctuations is characteristic of slow mag- 
netosonic modes, in which the cancellation is exact to 
~ 0{(3~^) when the plasma /3 ^ 1. When such modes 
are placed in a rotating shear flow whose rotation rate 
declines outward, they become the magnetorotational in- 
stability. We might naturally expect to see many such 
features in our shearing box simulations, and some can 
be expected to grow to nonlinear amplitude. 

Figure [TJ] shows correlation plots between magnetic 
pressure fluctuations, density fluctuations, and thermal 
(i.e. gas plus radiation) pressure fluctuations in the 
z = midplane of simulation 0519b near the 200 or- 
bit epoch. Most cells (i.e. most of the area as indicated 
by the red and yellow regions in all three panels) house 
weakly negative magnetic fluctuations, and the thermal 
pressure and density fluctuations in these cells approxi- 
mately obey the expected adiabatic relation for acoustic 
waves: (JPthorm/fthorm = Ti5p/p, where Fi ~ 4/3 is 
the flrst generalized adiabatic index for a gas and radi- 
ation mixture (Chandrasckhar 1967). However, in the 
upper right plot of magnetic pressure perturbation vs. 
thermal pressure perturbation, these acoustic waves os- 
cillate horizontally back and forth about a mean ther- 
mal pressure perturbation that is set by a background of 
isobaric fluctuations. This is indicated by the fact that 
the contours of cell-counts stretch diagonally across this 
plot, showing that at large amplitude the magnetic pres- 
sure and the total thermal (mostly radiation) pressure 
are anti- correlated. In addition, the magnitudes of the 
two fluctuations are similar, so that |(5Pmag + (^Pthcrml is 
in general considerably smaller than I^Pmagl + I'^Pthcrml- 
The isobaric modes also exhibit a clear anti-correlation 
between magnetic pressure and density, presumably be- 
cause the density is positively correlated with the gas 
pressure component of the thermal pressure. 

While most of the volume is occupied by weak mag- 
netic fleld. Figure [TSl shows that the regions of enhanced 
magnetic pressure associated with the isobaric modes are 
numerous enough to dominate the total magnetic energy 
budget in the vicinity of the midplane. As much as half 
the magnetic energy can be located in regions where the 
magnetic pressure is more than twice the local horizon- 
tal average. The maximum enhancement ratio in a single 
horizontal slice can be as high as ~ 10 at certain times 
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Fig. 11. — Spatial distribution of various quantities in the 2 = midpl ane of simulation 0519b at epoch 200.3 orbits (near the time of peak 
amplitude of the = 9, Uy = I nonaxisymmetric wave shown in Figure [T3t . On top from left to right are the total (radiation plus gas plus 
magnetie) pressure perturbation, scaled with the horizontally averaged total pressure at the midplane; the fluid velocity divergence, scaled 
with the horizontally averaged midplane sound speed divided by the fiducial scale height; and the density perturbation scaled with the 
horizontally averaged midplane density. On the bottom from left to right are the radiation pressure perturbation, gas pressure perturbation, 
and magnetic pressure perturbation, respectively, scaled with the horizontally averaged total pressure at the midplane. 



and places. 



5. RADIATIVE DAMPING 



Any phenomenon involving radiation pressure fluctu- 
ations (acoustic or isobaric) can be dissipated (i.e., lose 
its mechanical energy to heat, with associated entropy 
production) via photon diffusion, which enters our for- 
malism in the diffusion equation. In fact, we find that 
both types of radiation pressure fluctuation contribute 
significantly to the radiative damping, although the iso- 
baric perturbations dominate. 

One way of seeing this is to examine the integrated ra- 
diative damping as a function of fractional magnetic and 
thermal pressure fluctuations. This is shown in the up- 
per left plot of Figure [16] for the 200 to 202 orbit interval 
at the z = midplane in simulation 0519b. (The cor- 
responding incidence of fractional magnetic and thermal 
pressure fluctuations themselves was shown in the upper 
right hand plot of Figure [TH) Two distinct regions in 
this plane contribute to the radiative damping, one asso- 
ciated with large positive magnetic fluctuations stretch- 
ing up the isobaric locus, and one associated with slightly 
negative magnetic fluctuations that can occur with both 
positive and negative thermal pressure fluctuations that 
are mostly due to acoustic waves. 

Integrating the data of the 2-d distribution function 
over all thermal pressure perturbation values results in 
the distribution shown in the upper right plot of Fig- 



ure [TBI Because weak magnetic fields (negative magnetic 
pressure perturbations) dominate the spatial volume in 
the midplane, radiative damping from negative magnetic 
pressure perturbations here is primarily due to acoustic 
waves. In contrast, positive magnetic pressure pertur- 
bations are mostly associated with the strongly magne- 
tized, nonlinear isobaric fluctuations. The contributions 
to radiative damping from these two types of fluctua- 
tion are comparable: 35 percent for acoustic waves, and 
65 percent for isobaric fluctuations for the 200-202 orbit 
epoch in simulation 0519b, for a total midplane radiative 
damping rate of 9 x 10^'* ergs cm""^ s~^. 

Because the acoustic waves are much faster oscillat- 
ing fluctuations, they are presumably responsible for the 
rapid variability seen in the time dependence of the box- 
integrated radiative damping (upper dashed gray curve) 
of Figure [9] In fact, the origin of the high frequency 
spikes in the power spectrum of the box-integrated ra- 
diative damping shown in Figure [10] is now clear. All of 
these spikes are higher order harmonics of a fundamen- 
tal frequency of ~ 2.36 inverse orbital periods. This fre- 
quency is exactly the inverse of the shear time of the box, 
and reflects the time interval between successive pulses 
of the dominant |7ij,| — 1 nonaxisymmetric waves (see 
Figure [111). 

The measured contributions to the radiative damp- 
ing can also be estimated analytically from the observed 
fluctuation amplitudes, at least for the nonaxisymmetric 
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Fig. 12. — Spatial distribution of various quantities in tlie i; = vertical plane of simulation 051 9b at epoch 200.3 orbits, depicting only 
the regions within ±0.7H of the z = midplane. The quantities shown are the same as in Figure [TT] The perturbations are all measured 
with respect to, and scaled by, locally horizontally averaged quantities at each height z. 
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Fig. 13. — Time dependence of the real part of the Uy = 1 spatial 
Fourier amplitude in total pressure at the z = midplane, scaled 
with the horizontally averaged total pressure there. Different col- 
ored curves correspond to different values of n^;: = 5 (red), 
nx = 6 (green), rix = 7 (blue), = 8 (yellow), = 9 (cyan), 
and Tlx = 10 (magenta). Each of these values correspond to 
wave vectors that happen to be swinging from leading to trailing 
in or near the particular time interval shown. Other values of 
are plotted as black curves. 



fast waves, which are not too nonlinear. All the fluctua- 
tions that we observe in the midplane regions have length 
scales such that the photon diffusion time is much longer 
than the sound crossing time. When radiation pressure 
dominates both magnetic and gas pressure, the exponen- 
tial damping rate of the amplitude of a linear, plane wave 
fast mode (acoustic wave) in this slow diffusion limit is 



given by 



(23) 



where k is the wavenumber of the mode. Physically, 
acoustic waves in a radiation dominated plasma are sim- 
ply damped at the rate that photons diffuse across a 
wavelength. The acoustic radiative damping rate can 
therefore be estimated by multiplying the average wave 
energy density by twice this damping rate (twice because 
energy is proportional to the amplitude squared) , i.e. 
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(24) 



where ct ~ 6 x 10^ cm s ^ is the total sound speed in 



the radiation dominated plasma (see Figure . From 
Figure 1111 we estimate a typical pressure amplitude 
of <5Ptot,inax/^tot ^ 0.06 and a wavenumber of fc ^ 
27r/(0.3i?). For these numbers, we find Qrad, acoustic ~ 
3 X 10^** ergs cm~'^ s~^, very close to our measured 
acoustic contribution near this epoch of 35 percent of 
9 X 10^'' ergs cm~'^ s~^. 

The short two-orbit epoch we have analyzed in detail in 
this section is not atypical for this simulation: the mid- 
plane radiative damping rate of 9 x 10^^ ergs cm~^ at 
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Fig. 14. — Two dimensional distributions of the incidence of given values of magnetic pressure fluctuations, density fluctuations, and 
thermal pressure fluctuations in the 2 = midplane, averaged over the 200 to 202 orbit time period in simulation 0519b. Each distribution 
is normalized such that the two dimensional integral over the perturbation variables is unity. The red line in the upper right plot shows 
the expected relation for isobaric fluctuations, i.e. where the sum of the magnetic and thermal pressure perturbations vanish. The red line 
in the lower left plot shows the expected relation for adiabatic perturbations: (Pthorm ^ Ahcrm,av)/^thcrm,av = ^liP ~ PsLv)/pa.v 

o_sr^ — ' — ' — I — ' — ' — ' — I — ' — ' — ' — I — ' — ' — ' — I — ' — ' — ^ Although the contributions of acoustic waves and iso- 

baric fluctuations to radiative damping are comparable 
near the midplane, plots similar to Figure [12] at fixed 
heights more than ~ 0.5H away from the midplane 
clearly show that isobaric fluctuations become dominant 
at these higher altitudes. At these locations, the two- 
dimensional distribution of radiative damping as a func- 
tion of magnetic and thermal pressure closely follows the 
isobaric locus for both positive and negative magnetic 
perturbations, in contrast to the midplane distribution 
(upper left plot of Figure \W\\ . This fact is consistent 
with the patterns of fluctuations that we observe as a 
function of height in Figure 1121 In that flgure, we see 
that the vertical yellow and blue stripes near the mid- 
plane associated with nonaxisymmetric acoustic waves 
in total pressure (upper left plot) extend out only to 
\z\ ^ 0.2-0.3H. Beyond that, isobaric fluctuations of 
both signs of thermal and magnetic pressure dominate 
the volume (lower left and right plots). 




Fig. 15. — The differential distribution of fractional contributions 
to the total magnetic energy at a given height, as a function of 
local magnetic pressure scaled by the average magnetic pressure at 
that height. The distributions are averaged from 200 to 202 orbits 
in simulation 0519b, and are plotted for the midplane (2 = 0, 
black), 2 = —0.5H (red), and 2 = +0.5-ff (green). Each of the 
distributions is normalized such that the integral with respect to 
scaled magnetic pressure is unity. 



this time is close to the midplane radiative damping rate 
time-averaged over the entire ~ 400 orbits of the simula- 
tion (upper gray dashed curve of Figure El). We reiterate 
here how significant these numbers are (as much as 22 
percent of the work done by the walls when integrated 
over the entire box) for the overall energy budget at these 
high levels of radiation pressure support. 



6. THE NATURE OF VERTICAL ADVECTIVE FLUXES 

We now turn to examine two possible origins of the 
vertical radiation advection energy flux: transport by 
trapped vertical waves and buoyancy. 

6.1. Vertical Epicyclic and Acoustic Waves 

Figure [T7] shows the temporal power spectrum of the 
vertical radiation advection energy flux at every height z 
in simulation 0519b. A hierarchy of vertical modes is ap- 
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Fig. 16. — Upper left: two dimensional distribution function of the z = midplane radiative damping rate as a function of fractional 
thermal pressure perturbation and fractional magnetic pressure perturbation, averaged over the 200 to 202 orbit interval in simulation 
0519b. The distribution is normalized such that the integral over both perturbation variables is unity. The red line shows the locus of 
perfect isobaric perturbations. Upper right: projection of the distribution of radiative damping across all thermal pressure perturbations, 
as a function of magnetic pressure perturbation. The regions largely associated with acoustic waves and isobaric fluctuations are indicated. 
Lower left: projection of the distribution of radiative damping across all magnetic pressure perturbations, as a function of thermal pressure 
perturbation. 
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Fig. 17. — Vertically resolved power spectrum of the horizontally- 
averaged radiation advection flux (i?i'z)av in simulation 0519b. 

parent at frequencies at and above the orbital frequency. 
All but the lowest frequency of these modes are identical 
to those causing the oscillations and discrete frequency 
peaks in the mechanical pumping shown in Figures IHl and 

m 

The lowest order mode at the orbital frequency is just 
a vertical epicyclic oscillation, excited by vertical asym- 
metries in momentum losses through the vertical bound- 
aries of the box. The vertical velocity Vz in this mode is 
the same at all heights, and the mode simply uniformly 



displaces the entire plasma up and down in the box. It 
produces zero compression or rarefaction, and therefore 
does no pressure work, which explains its complete ab- 
sence in the temporal power spectra of pressure work 
terms shown in Figure [TUl Because it bodily displaces 
trapped photons up and down, it modulates the radia- 
tion advection energy flux, but it produces zero radiation 
advection through the plasma. 

The next lowest frequency mode is a vertical breath- 
ing mode, and higher and higher frequency modes clearly 
represent standing acoustic waves with increasing num- 
bers of vertical velocity nodes. The vertical velocities of 
these modes also modulate the radiation advection, but 
the net time-averaged flux would be zero if these modes 
were purely adiabatic. However, radiative damping of 
these modes does produce a net radiation advection en- 
ergy flux through the plasma over time. With a small 
amount of analytic work (see appendix for details), we 
can directly calculate this flux and compare it to the 
measured value. 

The strongest mode is the breathing mode, and we 
measured its velocity amplitude as follows. First we com- 
puted unwindowed, unbinned power spectra of the ver- 
tical velocity time series at each height in the box. We 
then summed the power at each height over the measured 
mode line profile in these power spectra (corresponding 
to the lowest peak in the upper green and black curves 
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of Figure [10)) , and fit the resulting power at each height 
in the midplane regions with the square of the expected 
mode eigenfunction (hnear for the breathing mode). Fi- 
nahy, we used Parseval's theorem to determine the ac- 
tual velocity amplitude of the mode. The results were 
~ (1 to 7) X 10® cm s~^ (roughly one percent of the ther- 
mal sound speed) across the simulations at one pressure 
scale height from the midplane. 
The damping rate of the breathing mode is (see eq. 

Trad -jr, (25) 

where Eg is the radiation energy density at the midplane. 
Not surprisingly, this is approximately the reciprocal of 
the nominal thermal time of the diskj^ Because the 
breathing mode is the longest vertical wavelength acous- 
tic mode, it should be damped at a rate given roughly 
by the rate at which photon diffusion can cool the disk. 

The radiative damping of the breathing mode causes 
phase changes between the velocity and pressure per- 
turbations, which result in a secular flux of mechanical 
energy carried by the mode. This energy flux is given by 
(see eqs. IA19l and lA20l in the appendix) 



Sv 



therm 



A I 5! 



"9~ 



All 
(26) 

where A\, is the velocity amplitude of the breathing mode 
at one pressure scale height Hp from the midplane. 

Evaluating equations (|25l) and (|26l). we find a peak ra- 
diation advection flux of only 7 x 10^® ergs cm"^ s~^ 
at the breathing mode amplitude observed in the 0519b 
simulation; the other simulations have comparable fluxes, 
although the ones with lower radiation to gas pressure ra- 
tios are smaller by factors of several. However, Figured] 
shows that the actual time-averaged peak radiation ad- 
vection flux in 0519b is more than five orders of magni- 
tude larger: 1 — 2 x 10^^ ergs cm"^ s~^. The breathing 
mode amplitude would have to be at least ~ 400 times 
larger, i.e. several times larger than the sound speed, 
in order to produce that large a flux. Similar discrepan- 
cies exist in all six simulations, and we therefore conclude 
that the acoustic waves do not contribute significantly to 
the time-averaged outward radiation advection. Instead, 
they simply modulate it at high frequencies. 

Similar calculations can be used to estimate the dis- 
sipation rate due to radiative damping of the breath- 
ing mode. For example, in simulation 0519b, this is 
~ 4 X 10^° ergs cm~^ s~^ at the midplane, orders of mag- 
nitude smaller than the actual rates of radiative damping 
and pressure work that we find in the midplane regions 
of this simulation (Figure [6]) . The same is true for all 
the other simulations, and we conclude that the stand- 
ing vertical acoustic waves contribute negligibly to radia- 
tive damping and the time-averaged mechanical pumping 
work, although they do significantly modulate the latter. 

^ The actual thermal time as measured by the time average 
of the box-averaged thermal energy content divided by the sum 
of the top and bottom emergent radiation fluxes is several times 
shorter than this nominal thermal time, presumably because the 
dissipation profiles peak off the midplane and some of the energy 
is carried outward by mecha nical motions (see discussion around 
eq. A8 of lHirose. Blaes. k. KroUk 200d ). 
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t (orbits) 

Fig. 18. — Itorizontally-averagcd radiation advection flux 
(£^fz)av, smoothed over a running two orbit interval, as a func- 
tion of height z and time t in simulation 0519b. Fine scale vertical 
striations in this plot are what remain of the vertical epicyclic and 
acoustic modulation of the radiation advection after the two-orbit 
smoothing. 
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Fig. 19. — Horizontally-averaged Poynting flux, smoothed over a 
running two orbit interval, as a function of height z and time t in 
simulation 0519b. 



6.2. Buoyancy as the driver of advection 

We have just seen that, while the radiation advec- 
tion is modulated by vertical acoustic modes, these 
waves do not contribute significantly to the net radi- 
ation advection flux. If we smooth over this rapid 
modulation, a clear pattern emerges as shown in Fig- 
ure 1181 The radiation advection flux is modulated on 
long (~ 5 orbit) time scales, and this modulation is 
in fact highly correlated with a similar modulation in 
the smoothed spacetime pattern of vertical Poynting flux 
shown in Figure 1191 These spacetime patterns of Poynt- 
ing flux have been seen in all previous vertically strat- 
iflcd MRI simulations, both with and wit h out explicit 
thermodynamics (.Brandciiburg et al.l 119951: IStone et al.l 



19961: IMiUer fc Stondl2000l:lHirose. Krolik. fc Stonel l2006r 
Davis et al.ll2009l: iShi, Krolik. fc Hirosdl2010[ ). but their 
origin has remained a mystery. 

We argue here that both the radiation advection and 
the vertical Poynting flux are due to buoyancy in the re- 
gion inside \z\ ^ H . That there is any buoyancy at all in 
the midplane regions is at first surprising, as the time and 
horizontally averaged pressure and density profiles are 
linearly stable to buoyant instabilities in the midplane 
regions. Indeed, when the square of the hydrodynamic 
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Brunt- Vaisala frequency, 
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dlnp 

dz 



(27) 



is computed using horizontally- and time-averaged pres- 
sure and density profiles, it is positive everywhere off the 
midplane in all our simulations, indicating a hydrody- 
namically stable average vertical entropy profile. 

Magnetic fields can still cause buoyancy instabilities 
in the form of the interchange and undulatory Parker 
modes. The interchange mode is stable provided the ra- 
tio of magnetic field strengt h to mass densit y does not de- 
crease outward too fast (e.g. IAchesonlll979[ ). In our simu- 
lations, the time-averaged profiles have d/dz \yv{B/ p) > 
everywhere off the midplane, so the interchange mode is 
linearly stable. 

This leaves the undulatory Parker modes, which 
are typically linearly unstab le in the surface layers 
(|Blaes. Hirose. fc Krolikir2007D . In the midplane regions, 
where radiative diffusion is slow, the linear stability crite- 
rion (jNewcomh. 19611 can be expressed in terms of the re- 
quirement that the square of the magnetic Brunt- Vaisala 
frequency be positive. This can be defined as 
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(28) 



where ct = (PiPthcrm/p)^^^ is the adiabatic sound speed 
in the gas plus radiation mixture, wa is the Alfven speed, 
and the last equality follows if the field is purely horizon- 
tal with no vertical tension forces. The solid lines in Fig- 
ure [201 show that the midplane regions are stable to the 
undulatory Parker modes by this criterion. If radiative 
diffusion is fast enough to suppress temperature fluctua- 
tions, the undulatory Parker stability criterion (Gilmani 
can instead be written in terms of a modified mag- 
netic Brunt- Vaisala frequency 



4) 



dz 



(29) 



where q = {p/pY^"^ is the isothermal sound speed in 
the gas. The dashed lines in Figure BUI show the profiles 
of N"^^^ .^. Hence even if radiative diffusion were fast 
in the midplane regions, they would still be stable sim- 
ply because the time and horizontally-averaged magnetic 
pressure typically increases outward in these regions. 

The midplane regions, viewed in terms of horizontally- 
averaged quantities^ are therefore linearly stable to buoy- 
ant perturbations. Nevertheless, finite amplitude pertur- 
bations with three dimensional structure within the tur- 
bulence can still be locally buoyant. As we have already 
seen, magnetosonic slow modes create regions of high 
magnetic pressure and low gas density. Here they have 
nonlinear amplitude, but retain these qualitative charac- 
teristics. Because they have the same pressure as other 
locations at that altitude, but have lower density, they 
are, of course, buoyant. Such low density, high magnetic 
pressure fluctuations are clearly present throughout the 
region within ±0.7iJ of the midplane shown in Figure [121 

Unlike the breathing modes, these fluctuations do carry 
enough energy flux to explain the advection. Figure 
shows the two dimensional distribution of radiation ad- 
vection flux at z = +0.5H in simulation 0519b, aver- 
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Fig. 20. — Square of magnetic Brunt- Vaisala frequency appro- 
priate for the undulatory Parker instability, as computed from the 
time averaged pressure and density profiles of simulations 1112a 
(top) and 0519b (bottom). T he s olid curves are for slow radiative 
diffusion (first equality of ea^ l28l l. while the dashed curves are for 
rapid radiative diffusion (eq. [29J. (The second equality of equation 
1281 produces curves that deviate from the solid curves only in the 
surface regions where magnetic tension forces are significant.) 

aged over 200 to 208 orbits, as a function of magnetic 
pressure and density perturbations. Integrating over all 
densities and magnetic fields to compute the total net 
radiation advection gives 9 x 10^^ ergs cm~^ s~^, which 
is of course consistent with the horizontally and time- 
averaged radiation advection at this height around this 
epoch (Figure fTS]). It is also consistent with the average 
radiation advection flux that we see at this height over 
the entire duration of the simulation (Figure [ij . 

Most importantly, it is clear that magnetic pressure 
and density continue to be anti-correlated above the mid- 
plane, just as they are in the midplane. Moreover, out- 
ward radiation advection tends to be associated with low 
density and high magnetic pressure, and inward radiation 
advection tends to be associated with high density and 
low magnetic pressure. This can also be seen if we in- 
tegrate over all densities (the distribution in the upper 
right plot) or all magnetic pressures (the distribution in 
the lower left plot). Thus, the advection is clearly asso- 
ciated with these local regions of buoyancy. 

A corollary of this flnding is that the regions respon- 
sible for the bulk of upward radiation advection actually 
have radiation energy density somewhat smaller than 
other regions at the same altitude. In other words, this 
process, although it bears some resemblance to classical 
convection, differs strikingly in that it carries heat up- 
ward in fluid elements that are initially cooler than their 
surroundings. Instead of being due to a higher entropy 
per unit mass (which in fact they have in spite of being 
cooler than their surroundings) , their buoyancy is driven 
by a higher magnetic pressure per unit mass. 

The fact that cooler than average buoyant fluid ele- 
ments can still give rise to a net outward energy flux can 
also be seen by the following simple argument (3 In any 
given horizontal plane, divide the radiation energy den- 
sity, mass density and vertical velocity into their horizon- 
tal average plus perturbations (just as we did with other 
variables to separate radiative damping from mechanical 
pumping in the pressure work terms in section lX^ above) : 
E = Sav + 6E, p = pav + Sp, and = Wz.av + Svz- We 

^ We thank the referee for suggesting this argument to us. 
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Fig. 21. — Upper left: two dimensional distribution of the 200 to 208 orbits time averaged radiation advection at 2 = +0.5H in simulation 
0519b, as a function of scaled magnetic pressure perturbation and density perturbation. The distribution is normalized such that the two 
dimensional integral over both perturbation variables is unity. Upper right: projection of the radiation advection distribution across all 
density perturbations, as a function of magnetic pressure perturbation. Lower left: projection of the radiation advection distribution across 
all magnetic pressure perturbations, as a function of density perturbation. 



observe zero net vertical mass flux in the simulations, so 
Vz,a,v = —{SpSvz)a,v/piiv The horizontal average of the 
radiation advection energy flux is therefore given by 



6Vz 



6E 



Pav 



(30) 



The fractional density fluctuations are generally much 
larger in magnitude than the fractional temperature (or 
radiation pressure or energy density) fluctuations, as we 
can see in Figures [TT] and [l2l It is the sign of these 
density fluctuations that give rise to an overall outward 
mean radiation advection energy flux. Put another way, 
the underdense regions give rise to a net outward horizon- 
tal average velocity, and because the temperature fluctu- 
ations are so small, this net outward velocity times the 
average radiation energy density is approximately the net 
outward radiation advection flux. 

While the rising magnetized fluid parcels start out 
cooler than their surroundings, they nevertheless even- 
tually become hotter than their surroundings, and it is 
for this reason that they are able to transfer heat from 
the hot midplane regions to the cooler regions away from 
the midplane. We can see this qualitatively through 
the following linearized treatment. Consider a bundle 
of (mostly horizontal) fleld lines that rises an inflnitesi- 
mal height Az from some initial height z. We take the 
gas pressure to be negligible, and assume that radiative 
diffusion is slow so that the plasma inside the bundle un- 



dergoes an approximately adiabatic change (as shown by 
Fig. [231 in this region the diffusion speed on a lengthscale 
0.1-ff is about comparable to the typical upward advec- 
tive speed in this region). We further assume that the 
mean gradient in the total pressure (radiation plus mag- 
netic) is in hydrostatic balance with gravity. Pressure 
equilibrium with the surroundings then implies that the 
radiation pressure P*^^ of the plasma inside the bundle 
will change by 
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(31) 



where P^^g is the initial magnetic pressure inside the 
bundle. On the other hand, the radiation pressure in the 
local background will change by 
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(32) 



where dPmag/c^Prad is a derivative following background 
pressures as they vary with height. We therefore find 
that the rising fluid parcel will cool less rapidly with 
height than the background provided 
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(33) 



This is nearly always true in our simulations, as the bun- 
dle starts out more magnetized and cooler than its sur- 
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roundings, and the background magnetic pressure gen- 
erally falls less rapidly with he ight than the radiation 
pressure (see e.g. Fig. 16 of iHirose. Krolik. fc BlaesI 
[2009). In fact, near the midplane, the magnetic pres- 
sure generally rises with height, while the radiation pres- 
sure falls. For example, in the advective region of 0519b 



d In P„ 



— rises outward from 



-1.5 very 



d In Prad 

near the midplane to a plateau at ~ 0.5. Thus, rising 
magnetized fluid parcels eventually become hotter than 
their surroundings, at which point heat will be trans- 
ferred outward to the surroundings by radiative diffusion. 

Not surprisingly, regions of strong Poynting flux are 
even more strongly associated with high magnetic field, 
low mass density isobaric fluctuations, as shown in Fig- 
ure [211 This flgure is constructed in a fashion exactly 
analogous to Figure [51] illustrating radiation advection. 
Supporting our argument that the two are fundamen- 
tally the same, the patterns seen in this figure are very 
similar to those seen in the radiation advection flgure. 
The principal contrast is that, unlike the radiation ad- 
vection case, there is very little downward Poynting flux; 
it is almost exclusively upward. This follows, of course, 
from the fact that the upward-moving isobaric fiuctua- 
tions have greater than average magnetic energy density, 
whereas they have smaller than average radiation energy 
density. 

6.3. Characteristic speed of advective motions 

Because the outward advection fluxes are associated 
with regions of enhanced magnetic fleld, we can define 
a horizontally-averaged vertical advection speed by us- 
ing the z-component of the Poynting vector S = c(E x 
B)/(47r), 



vLy/2 



(34) 



We compare the time-averaged vertical profile of this 
speed in simulation 0519b to the time-averaged vertical 
profiles of the thermal sound speed Ct and the Alfven 
speed in Figure [23FI The advection speed is several 
percent of the sound speed in the near-midplane region 
where advection is most important, but rises outward, 
reaching Mach numbers a little less than unity in the 
outer layers of the disk. This outward increase in Wadv 
is consistent with the characteristic buoyant rise time of 
5-10 orbits for magnetic features: because H ^ ct/il, 
such a characteristic time implies a mean upward Mach 
number ^ 0.1-0.2. 

Figure[23|also shows the characteristic one-dimensional 
speed of radiation diffusion Vdm over a length scale of ^ = 

* This figure shows two other interesting features that are out- 
side the scope of the current discussion. First, the thermal sound 
speed is remarkably constant with height, a property that is shared 
by all the simulations we analyze in this paper. Second, the Alfven 
speed is not very much larger than the thermal sound speed in the 
surface layers in simulation 0519b. This is in marked contrast to 
all the other simulations that have lower radiation to gas pressure 
ratio, which have very stron gly magnetized surface regions that are 
Parker unstable HB laes, Hi rosc. & Krolik 2007). Parker instability 
dynamics is still apparent in the surface regions of 0519b, but it 
is conceivable that at still higher levels of radiation support that 
magnetic support and Parker dynamics would become less impor- 
tant. 



O.lff, the characteristic size of the buoyant fiuctuations 
(Figure [12]) . We define this speed as this length scale 
divided by the radiation diffusion time, i.e. 



(35) 



The diffusion speed increases rapidly away from the mid- 
plane as the density drops. It is comparable to the ad- 
vection speed in the midplane regions, and then greatly 
exceeds the advection speed at heights greater than ~ H . 
This coincides with the heights above which radiative dif- 
fusion finally becomes dominant over radiation advection 
in this simulation (bottom panel of Figure [T]) . 

The average advection velocity profile also allows us 
to estimate the work done by the plasma in driving ra- 
diation advection. Figure [23| shows that at z = H , 
Uadv — 10^ cm in simulation 0519b. The average ra- 
diation pressure at that height is ~ 8 x 10^'' dyne cm~^. 
Hence the mechanical pumping work associated with ra- 
diation advection at this height is ~ —{E/'i)dva,dv/dz ~ 
-{E/3)va_dv/H ^ -2 X 10^^ erg cni^^ s^^. This agrees 
very well with the actual mechanical pumping in Figure[6] 
(lower gray dashed curve) . We conclude that the vertical 
expansion motions associated with radiation advection 
are largely responsible for the mechanical pumping work 
done by the plasma. 

7. DISCUSSION 

We have now answered a number of the questions posed 
at the beginning of this paper on the thermodynamics 
of radiation pressure dominated disks. Central to our 
discussion is the existence of high magnetic energy den- 
sity, low gas density regions that manifest within the 
turbulence as nonlinear slow magnetosonic modes. We 
named these regions "isobaric fluctuations" because the 
total pressure perturbation within them (thermal plus 
magnetic) is quite small. Higher than average magnetic 
pressure is associated with lower than average thermal 
pressure. The isobaric fluctuations are nonlinear, in that 
they are associated with order-unity magnetic intermit- 
tency. The magnetic energy density in these fluctuations 
is a few times the volume mean, and they constitute a 
signiflcant fraction of the magnetic energy in the turbu- 
lence. On the other hand, these regions are rarely of 
sufficient size and coherence that they can be thought of 
as "flux tubes" . 

We have elucidated the nature and role of radiative 
damping. To supply the energy that is ultimately dis- 
sipated by photon diffusion, thermal pressure does work 
on the fluid in association with compressive fluctuations. 
This can be distinguished from work done by the fluid 
in driving vertical mechanical motions by subtracting off 
appropriate horizontal averages in presure and velocity 
(equation [20]) . The compressive fluctuations responsi- 
ble for radiative damping are a combination of isobaric 
fluctuations and acoustic waves. In the midplane re- 
gions of our most radiation dominated simulation, we 
find that the dissipation associated with radiative damp- 
ing is roughly 65% due to isobaric modes and 35% to 
acoustic waves. At higher altitude, isobaric fiuctuations 
dominate the dissipation associated with radiative damp- 
ing- ^ 

It is noteworthy that the time-averaged spatial distri- 
bution of radiative damping is similar in shape to the 
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Fig. 22. — Upper left; two dimensional distribution of the 200 to 208 orbits time averaged Poynting flux at z = +0.5H in simulation 
0519b, as a function of scaled magnetic pressure perturbation and density perturbation. The distribution is normalized such that the two 
dimensional integral over both perturbation variables is unity. Upper right: projection of the Poynting flux distribution across all density 
perturbations, as a function of magnetic pressure perturbation. Lower left: projection of the Poynting flux distribution across all magnetic 
pressure perturbations, as a function of density perturbation. 
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Fig. 23. — Time- and horizontally-averaged vertical profiles of 
thermal sound speed (black), Alfven speed (solid blue), radia- 
tion diffusion speed over a length scale of 0.1_ff (solid red), and 
magnetically- weighted upward advective speed (dashed blue), for 
simulation 0519b. 

time-averaged spatial distribution of magnetic and ki- 
netic energy dissipation (Figures [5] and [51) ■ Similarly, the 
box-averaged radiative damping is well correlated with 
the magnetic and kinetic energy dissipation rates as a 
function of time: averaging over the short time scale fluc- 
tuations associated with acoustic waves in the radiative 
damping (upper dashed gray) curve of Figure [HI makes 
this history's relation to the magnetic (blue) and kinetic 
energy (green) dissipation curves clear. In addition, as 
we have already noted, a significant fraction of the total 



magnetic energy is found in regions where the magnetic 
energy density is of order twice the mean, places where 
strong isobaric fluctuations have led to field concentra- 
tion. We therefore suspect that isobaric fluctuations play 
a significant role in the grid scale magnetic and kinetic 
energy dissipation. In the magnetic case, this is perhaps 
not surprising, as the large concentrations of magnetic 
field must be associated with sharp gradients of magnetic 
field, as well as high electric current densities. 

Surprisingly, we have also seen that when the ratio of 
radiation pressure to gas pressure is large enough, the to- 
tal dissipation rate due to radiative damping can become 
comparable to the gridscale loss of magnetic energy. If 
this trend continues to still larger radiation pressure, it 
will mean that we will have the unique advantage of be- 
ing able to resolve the principal scale of dissipation in 
astrophysical MHD turbulence. That will finesse a great 
many of the uncertainties and difficulties that stand in 
the way of a clear interpretation of MHD turbulence in 
many other contexts. 

Indeed, there is growing evidence that the saturated 
state of MRI turbulence may depend on the microscopic 
visco s ity and resistivity of the plasma dF romang et alj 
2007t ILesur fc Longarettil 120071: Isimon fc Hawlcv 200^ 
Davis et al.l I2009j r In particular, much attention has 
been paid to the sustainability and properties of the 
turbulence for different magnetic Prandtl numbers (the 
dimensionless ratio of viscosity to resistivity), because 
this sets the relative sizes of the velocity and magnetic 
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dissipation scales. The inner, radiation dominated re- 
gions of optically thick, geometrically thin accretion disks 
around black holes generally have magnetic Prandtl num- 
ber greater than unity if the viscosit y and resistivity 
are d ominated by Coulomb scattering (|Balbus fc Henril 
120081 ). In other words, the viscous scale, where veloc- 
ity fluctuations are damped, exceeds the resistive scale 
where magnetic fluctuations are damped. Insofar as ra- 
diative damping also suppresses velocity fluctuations, al- 
beit only those which are compressive, it reinforces this 
ordering: velocity fluctuations are damped on larger spa- 
tial scales than magnetic fluctuations in radiation domi- 
nated accretion disks. 

Perhaps the most important question that we posed at 
the beginning of this paper was how a radiation dom- 
inated disk handles a dissipation rate that exceeds the 
critical value of Q* = cfl^/K associated with hydrostatic 
and radiative equilibrium. As shown in Figures [5] and [SI 
such an excess always occurs in the midplane regions. We 
find that the radiative diffusion flux is indeed constrained 
by hydrostatic equlibrium to have a divergence equal to 
Q* , and excess energy flux is carried outward by radia- 
tive advection, fluid motion carrying trapped photons. A 
smaller (by energy flux is conveyed electromag- 

netically by the same fluid motions. At higher altitudes, 
the vertical energy flux becomes almost entirely diffusive 
radiation flux, simply because it becomes progressively 
harder for small fluid elements to trap photons as the 
density falls. 

The existence of substantial vertical energy flux carried 
by something other than radiative diffusion or Foynting 
flux is a signiflcant addition to the repertory of accre- 
tion disk physics. This is the first time that advected 
radiation flux has been shown to be important in disks. 

We emphasize that these advected energy fluxes are not 
classical convection, as the horizontal and time-averaged 
entropy profiles are stably stratified. This is also true 
of the horizontal and time-averaged magnetic profiles, at 
least in the midplane regions. These advected fiuxes are 
also not in any way related to the old argument that ra- 
diation dominated disks have constant density and are 
therefore convectively unstable. As we discussed in the 
introduction, this argument rested on the assumption 
that the dissipation per unit mass is constant. This is 
simply not true, nor is the density constant with height. 

Local buoyancy is the mechanism responsible for driv- 
ing both radiation advection and upward Foynting flux, 
and this buoyancy arises from the low densities of the 
nonlinear, highly magnetized isobaric fluctuations. Be- 
cause the lower than average thermal pressure within 
these fluctuations is mostly radiation pressure, these 
buoyant fluid parcels are cooler than their surroundings, 
at least initially. Nevertheless, they advect radiation in- 
ternal energy outward because their magnetic pressure 
makes them buoyant. The work done by the plasma in 
driving this local vertical expansion explains the mechan- 
ical pumping work that we separated from the radiative 
damping in the total PdV work. 

Although these buoyant regions play an important role 
in the energy budget in very radiation-dominated disks, 
they are also important for a different reason in disks 
with weaker radiation content. The same isobaric fluc- 
tuations can also explain the characteristic upwelling of 
magnetic flux seen in all stratified disk simulations. Ra- 



diation forces per se are unnecessary: the only prereq- 
uisite is that there be thermal pressure fluctuations ac- 
companied by proportionate density fluctuations that are 
balanced by opposite sign fluctuations in the magnetic 
pressure. Data from our earlier gas-dominated simula- 
tions show essentially the same dynamics as the data 
from radiation-dominated simulations that has played 
the primary role in this paper's analysis. 

The quasi-periodic build up of magnetic flux in the 
midplane regions is also associated with azimuthal fleld 
reversals, clearly indicating some sort of dyna.nio activity 
(Brandenburg ct al. 1995; Johanse n. Youdin. fc KlahiH 
|2009: Davis et al..,20q9^,Gressel,2010-iO'Nein et al.ll2010l: 
IShi. Krolik. fc Hirosell2010D . The fact that the resulting 
upwelling motions are energetically signiflcant is the flrst 
example that we know of in astrophysics of a magnetic 
dynamo that is important for the global energetics of the 
medium. For example, the solar dynamo is in part driven 
by the outward energy transport associated with convec- 
tion in the sun. In the case of a radiation-dominated 
accretion disk, the dynamo creates buoyant high mag- 
netic field regions which then play a significant role in 
energy transport by advecting radiation (and magnetic 
field) outward. We have seen that this radiation advec- 
tion is just as important as radiative diffusion in the mid- 
plane regions for the most radiation dominated simula- 
tion we have conducted so far. 

A final outstanding question is how the turbulence 
"knows" that it needs to produce a heating rate that is 
of order the critical value Q* required by hydrostatic and 
radiative equilibrium. We have seen that, in fact, when 
the overall radiation to gas pressure ratio is high, the 
dissipation rate generally exceeds Q* in the midplane re- 
gions, and the excess heat is carried outward by radiation 
advection, not radiative diffusion. But is this radiation 
advection therefore somehow regulated to match the re- 
quirements of overall thermal equilibrium? We suspect 
that the answer is no, at least not directly. 

The nonlinear isobaric fluctuations are an inherent fea- 
ture of intermittency in MRI turbulence and are related 
(somehow) to the dynamo process producing quasiperi- 
odic reversals in the azimuthal field. They are inher- 
ently buoyant, and therefore the resulting rate of radia- 
tion advection is simply proportional to how much radi- 
ation energy density is present in the fluid to be carried 
upward with the magnetic field. As we just discussed, 
these fluctuations do appear to be associated with grid 
scale dissipation that helps determine the local radiation 
energy density, but that is a slow process. Dissipation, 
which is still mostly magnetic in nature, has to occur over 
time scales comparable to the thermal time to build up 
enough heat to significantly change the radiation energy 
density, simply because the volume-averaged magnetic 
energy density is so small0 

Instead, any small excess radiation energy density that 
appears because the average dissipation exceeds Q* re- 
sults in an excess vertical diffusive radiation flux. The 
resulting excess radiation pressure force quickly (on the 
fast sound crossing time scale) produces a vertical ex- 
pansion of the medium. As a consequence, work is done 

^ This inherent time lag is ultimately w hy the disk is thermally 
stabl e in the radiation dominated regime ijHirose. Krolik. & Blaej 
l200gh . 
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by the fluid and this lowers the radiation energy density, 
restoring hydrostatic equilibrium. This mechanism may 
in part be responsible for the excitation of the standing 
vertical acoustic waves. It may very well be that hydro- 
static equilibrium cannot be so maintained if the turbu- 
lent dissipation becomes significantly stronger. Outflows 
may appear at levels of radiation dominance we have not 
yet explored in our simulations. 

Ultimately, the mechanism that sets the midplane dis- 
sipation to be in excess of Q*, and the level of ambient 
radiation energy density that is then advected outward 
by buoyant magnetic field, is the level of turbulence in 
the medium. However, what sets the saturation level of 
MRI turbulence is still an open question. 

Finally, we remark that photon bubbles have been pro- 
posed as an additional piece of physics that may al- 
ter t he vertical tran sport of energy at high luminosi- 
ties (jBegelmanl 120061 ). The characteristic length scale 
of these instabilities is ~ c?/(7, where q is the isothermal 
gas sound speed and g = is the acceleration due to 

gra vity ( Gamm ic 1998; Begclman 2001; Blacs & Socrates 
I2OO3I : [Turner et al.ll2005D . In simulation 0519b, this is 
less than the size of our grid zones at distances more 
than ~ Q.5H away from the midplane. Clearly these 
instabilities cannot be resolved in our simulations, and 
their possible impact on the thermodynamics remains a 
topic for future study. 

8. CONCLUSIONS 

The thermodynamics of a radiation dominated disk dif- 
fers significantly from that envisaged in standard static 
models of accretion disks. In addition to microscopic 
dissipation associated with the MRI turbulent cascade, 
radiative damping of compressive fluctuations plays an 
increasingly significant role as the radiation to gas pres- 
sure ratio increases. The compressive fluctuations asso- 
ciated with this damping consist of fast magnetosonic 
waves as well as nonlinear isobaric fluctuations arising 
from the turbulence itself. This damping is energetically 



significant: on the order of tens of percent of the over- 
all dissipation at the highest levels of radiation pressure 
support that we have simulated. It is also numerically re- 
solvable, in contrast to microscopic resistive and viscous 
dissipation. 

Buoyancy of the highly magnetized, low density iso- 
baric fluctuations is responsible in part for the butter- 
fly diagram that is always seen in simulations of MRI 
turbulence with vertical gravity. This buoyancy is in- 
trinsically three-dimensional in nature; even though the 
horizontally-averaged structure is very stably stratifled 
in the midplane regions, localized concentrations of mag- 
netic held generated by the turbulent dynamo are under- 
dense and still produce outward advection of magnetic 
held. When the plasma is radiation-dominated, these 
buoyant motions become significant for the overall en- 
ergetics of the plasma. Outward advection of photons 
becomes comparable to radiative diffusion, and the asso- 
ciated vertical expansion work must be included in bal- 
ancing vertical heat transport and dissipation. 
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APPENDIX 



VERTICAL ENERGY TRANSPORT BY STANDING VERTICAL MODES 

In this appendix we outline the derivations of the damping times and vertical energy fiuxes of the standing vertical 
acoustic waves discussed in section [6. II 

We consider a vertically stratified equilibrium with a background shear flow Vy — ~3Qx/2. Apart from Vy, all other 
fluid quantities are constant in the horizontal direction. We assume that the radiation and gas exchange heat sufficiently 
rapidly that Tj-ad — T everywhere. We also take the medium to be sufficiently optically thick that the radiation 
pressure tensor P only has diagonal elements i?/3, and radiation transport is diffusive with kP- ~ Kcs — constant. 
These thermodynamic assumptions are excellent approximations in the regions of the simulations where the acoustic 
waves are evident. Finally, we approximate the equilibrium magnetic field as being purely azimuthal and having only 
vertical gradients. 

The equilibrium vertical structure is given by hydrostatic equilibrium, 

radiative equilibrium, 

dz ^ 

and diffusive radiation transport. 



(A2) 



c dE 
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Here Pthcim = p+E/3 is the total thermal pressure. Note that we are including turbulence only insofar as it contributes 
energy dissipation Q. 

We now consider vertically propagating, longitudinal waves on this equilibrium. The Eulerian perturbed fluid 
velocity Sv{z,t) = d£^/dt, where ^ is the vertical Lagrangian displacement of each fluid element. The linearized 
continuity equation is 

where here and elsewhere, A = 6 + ^{d/dz) refers to Lagrangian perturbations. The linearized vertical momentum 
equation is 

The artificial viscosity term has been lost as it is nonlinear in the perturbation velocity (jStone fc Normanlll992[ ). The 
linearized flux-freezing equation is 

SB = -^{30 (A6) 
oz 

Finally, the linearized total internal energy equation may be written as 

I (i£f ) ~ ^4 {t) ^ (^^ " ^) ' ^^^^ 

where Fi an d F^ a re the usual generalized adiabatic exponents for a mixture of gas and radiation at the same 
temperature (jChand rasckhar f967), both of which are close to 4/3 for our radiation-dominated simulations. This 
equation may be integrated in time to give 

-fthorm P Ahcrm J \ OZ J 

Assuming the mode periods are short enough that there is little energy exchange between the modes and the 
turbulence, little turbulent dissipative heating, and little radiative damping, we can for the moment set 5Q = bF = 
in the equations. Assuming a time dependence oc exp(— zwi), equations (|A4[) - (|A6|) and (|A8[) may then be combined to 
give a single equation in Sv: 

d_ 

dz 



i l^thcrm + -; — -T^ — 

An I oz 



+ - = 0. (A9) 



The vertical epicyclic mode is the simplest general solution to this equation, with uj — and 5v being spatially 
constant. The breathing mode is the next simplest. Provided Fi is spatially constant (which is true if gas or radi- 
ation dominates the thermal pressure), and the magnetic energy density is much less than the thermal pressure (an 
assumption that typically fails only in the low density surface layers), then the hydrostatic equilibrium equation (jAl|) 
guarantees that w = (1 -I- TiY^'^Vl and 5v o^ z satisfy equation (IA9|) . 

Like the breathing mode, higher order modes are acoustic in nature , but their frequencies and eigenfunctions depend 
more sens itively on the details of the equilibrium vertical structure (jOkazaki. Kato k, Fukuelll987l : iLubow fc Pringld 
[1993; Blae s. Arras, fc Fr agile 2001). 

We now consider radiative damping of these modes. Multiplying equation (jA5l) by 5v, and then combining with 
equations (|A4p . (|A6[) . and (|A8[) . we finally obtain after some algebra the following wave energy conservation equation 



d_ 
di 



d I f B 

h -77- (5vAPthcrm + 5vA\ — 

OZ I \87r 
(F3 - 1) dAp /■* / d6F\ 



p dt 



The right hand side of equation (|A10[) represents sources and sinks of wave energy due to perturbed turbulent 
dissipation and radiative losses. True microscopic viscosity and resistivity would also contribute additional terms on 
the right hand side arising from additional terms in the perturbed momentum and flux-freezing equations, but our 
numerical dissipation scheme has no such terms. In reality they would be present, but likely even more important 
than these would be terms representing turbulent fluctuations, including turbulent viscosity and resistivity, which we 
have ignored here. However, we believe the damping to be dominated by radiative diffusion. 

As suming that the damping rate is small, we can estimate it by computing a work integral for the modes (e.g. iCoxl 
Il980f ). Approximating all the perturbed quantities on the right hand side as being perfectly periodic at the mode 
period H = 27r/wR, where wr, = Re(a;), we can time-average equation (jAlOl) over this mode period and integrate the 
right hand side by parts. Then integrating over all height, assuming negligible wave energy fluxes leaving the top and 
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bottom boundaries, we finally obtain for the exponential damping rate F = Im(a;) of the mode amplitude 



where 



^2 



(A12) 



is the mode energy. 

Because Sv and therefore ^ are spatially constant for the epicyclic mode, equation (jA4l) implies that Ap — 0, i.e. 

fluid elements preserve their density under this mode. The work integral therefore vanishes, and there is no damping 

of this mode. This is as one would expect: the only way to change this mode is to exert an external force on the 

system, or eject momentum through the vertical boundaries. 
The breathing mode and higher order acoustic modes are damped, however. The linearized radiative diffusion 

equation is 

SF^-f'J^^-^'-^. (A13) 

If we consider a radiation dominated equilibrium, and neglect gas and magnetic pressure contributions to the mode 
dynamics, then equations (jASp and (jA13l) give a simple relationship between the perturbed radiative flux and the 
perturbed velocity. 

Setting 5v — Ah{z/Hp) coswr^ for the breathing mode, where Ah is the velocity amplitude at one pressure scale height 
Hp and wr = (7/3)^/^17, we find a radiative damping rate of 

r.ad - . „ p / dz, (A15) 

where Eq is the midplane radiation energy density, and the integral goes over the entire equilibrium, or at least that 
portion which is optically thick. Clearly this integral can be set equal to the pressure scale height Hp times some 
numerical constant. For example, for an n = 3 isentropic polytrope, we obtain 

315 / cn'^ \ cn^ ,^ 

Trad = ^. (A16) 

256 \HcsEq J HcsEq 

The existence of this damping produces nonzero time-averages in the wave energy fluxes because it changes the 
phase difference between the velocity perturbation 5v and the pressure perturbations APthcrm and A(i?^/87r). In the 
absence of an exact solution to the non-adiabatic eigenfunctions, we compute these pressure perturbations using the 
adiabatic equations. Setting 5v = Ah{z/ Hp)e~^* coswr^, then the instantaneous energy fluxes are 

(JuAPthcrm ^ ^v—- = p7Te~^^* coscjRtsinwRt cos^WRt + O ^ (A17) 

3 9ionHp \ ujB. J \^rJ 

and 

SvA — = cos ujRt sin ujRt cos^wr^ ] + 0{ ^] . (A18) 

\8n J AnujRHp \ wr / V'^r/ 

The first term in each of these expressions periodically changes sign, and a time-average of this term over one mode 
period results in positive or negative values depending on the phase of the time-interval chosen. This term therefore 
produces no secular energy flux. The second term does, and the resulting time-averaged fluxes are, to lowest order in 

r/wR, 

. . 'i.EzTAl , 

< (JwAPthorm >= XTT^TTTT (A19) 



and 



2m^Hl 



P2\\ ^.B'^zTAl 



( — ) ) ^ (A20) 



'p 



Note that these two equations are completely independent of whether the decay rate F was due to radiation damping 
or some other process. 
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